﻿// This file is part of GEOTRON.

// GEOTRON is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.

// The GEOTRON library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.

// You should have received a copy of the GNU General Public License
// along with GEOTRON.  If not, see <http://www.gnu.org/licenses/>.

// The Ordinance Survey include a PDF file on their website which 
// describes the process of conversion to and from their grid system.
// The algorithms described are implemented here.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Geotron
{
    /// <summary>
    /// Represents a UK grid reference
    /// </summary>
    public class UKGridReference
    {
        public double Northings;
        public double Eastings;

        double f_0 = 0.9996012717;
        double phi_0 = 49;
        double lambda_0 = 2;
        double e_0 = 400000;
        double n_0 = -100000;
        double a = 6377563.396;
        double b = 6356256.910;
        double rho;
        double nu;
        double eta2;
        double e2;

        double vii;
        double viii;
        double ix;
        double x;
        double xi;
        double xii;
        double xiia;

        /// <summary>
        /// Instantiates initial values of phi_0 and lambda_0
        /// </summary>
        private void PerformInPlaceConversions()
        {
            phi_0 = ((double)49 / (double)180) * Math.PI;
            lambda_0 = ((double)-2 / (double)180) * Math.PI;
        }

        /// <summary>
        /// Internal function implementing formula C3
        /// </summary>
        private void C2(double phi)
        {
            e2 = ((a * a) - (b * b)) / (a * a);
            nu = a * f_0 * Math.Pow((1 - (e2 * Math.Pow(Math.Sin(phi), 2))), -0.5);
            rho = a * f_0 * (1 - e2) * Math.Pow((1 - (e2 * Math.Pow(Math.Sin(phi), 2))), -1.5);
            eta2 = (nu / rho) - 1;
        }

        /// <summary>
        /// Internal function implementing formula C3
        /// </summary>
        private double C3(double phi)
        {
            double n = (a - b) / (a + b);
            return b * f_0 * (((1 + n + (1.2 * Math.Pow(n, 2)) + (1.2 * Math.Pow(n, 3))) * (phi - phi_0)) - (((3 * n) + (3 * n * n) + ((21 / 8) * n * n * n)) * Math.Sin(phi - phi_0) * Math.Cos(phi + phi_0))
                + ((((15 / 8) * n * n) + ((15 / 8) * n * n * n)) * Math.Sin(2 * (phi - phi_0)) * Math.Cos(2 * (phi + phi_0))) - ((35 / 24) * n * n * n * Math.Sin(3 * (phi - phi_0)) * Math.Cos(3 * (phi + phi_0))));
        }

        /// <summary>
        /// Internal function implementing formula C6
        /// </summary>
        private double C6()
        {
            return ((Northings - n_0) / (a * f_0)) + phi_0;
        }

        /// <summary>
        /// Retrieve the co-ordinates represented by this grid reference.
        /// </summary>
        public Location Coordinates
        {
            get
            {
                // Properly initialize all our stuff
                PerformInPlaceConversions();
                double phi_dash = C6();
                double m = C3(phi_dash);
                while ((Northings - n_0 - m) >= (0.01 * Math.Pow(10, -3)))
                {
                    phi_dash = ((Northings - n_0 - m) / (a * f_0)) + phi_dash;
                    m = C3(phi_dash);
                }
                C2(phi_dash);
                vii = Math.Tan(phi_dash) / (2 * rho * nu);
                viii = (Math.Tan(phi_dash) / (24 * rho * nu * nu * nu)) * (5 + (3 * Math.Pow(Math.Tan(phi_dash), 2)) + eta2 - (9 * (Math.Pow(Math.Tan(phi_dash), 2)) * eta2));
                ix = (Math.Tan(phi_dash) / (720 * rho * Math.Pow(nu, 5) ) ) * (61 + (90 * Math.Pow(Math.Tan(phi_dash), 2)) + (45 * Math.Pow(Math.Tan(phi_dash), 4)));
                x = 1/(Math.Cos(phi_dash) * nu);
                xi = Math.Pow((Math.Cos(phi_dash) * 6 * Math.Pow(nu, 3)), -1) * ((nu / rho) + (2 * Math.Pow(Math.Tan(phi_dash), 2)));
                xii = Math.Pow((Math.Cos(phi_dash) * 120 * Math.Pow(nu, 5)), -1) * (5 + (28 * Math.Pow(Math.Tan(phi_dash), 2)) + (24 * Math.Pow(Math.Tan(phi_dash), 4)));
                xiia = Math.Pow((Math.Cos(phi_dash) * 5040 * Math.Pow(nu, 7)), -1) * (61 + (622 * Math.Pow(Math.Tan(phi_dash), 2)) + (1320 * Math.Pow(Math.Tan(phi_dash), 4)) + (720 * Math.Pow(Math.Tan(phi_dash), 6)));

                var location = new Location();
                location.Mode = AngleType.Radians;
                location.Latitude = (phi_dash - (vii * Math.Pow((Eastings - e_0), 2)) + (viii * Math.Pow((Eastings - e_0), 4)) - (ix * Math.Pow((Eastings - e_0), 6)));
                location.Longitude = (lambda_0 + (x * (Eastings - e_0)) - (xi * Math.Pow((Eastings - e_0), 3)) + (xii * Math.Pow((Eastings - e_0), 5)) - (xiia * Math.Pow((Eastings - e_0), 7)));
                return location;
            }
        }

    }
}
